% Canny Edge Detection Implementation
% Author:   TuanND
% Date:     02/2013
% Note: method: {'interpolation', 'average', 'fast'}
% interpolation: 8 directions + interpolation for sub-pixel
% average: 8 directions + average for sub-pixel
% fast: 4 directions + neareast neigbour

function [canny_edge_img color_edge_img mg] = ecen5283_edge_detect_canny(src_img, sigma, hsize, method)
[rows cols] = size(src_img);
src_img = double(src_img);

%% generate first order differentation of Gaussian kernel
    function dgx = gen_dgauss_x(sigma, haftsize)
        dgx_size = 2 * haftsize + 1;
        dgx = zeros(dgx_size, dgx_size);
        for m = 1:dgx_size
            for n = 1:dgx_size
                u = m - haftsize -1;
                v = n - haftsize - 1;
                dgx(m,n) = (-u/(2*sigma^4)) * exp(-(u^2 + v^2)/(2*sigma^2));                
            end
        end                
    end
hx = gen_dgauss_x(sigma, hsize);
hy = hx';

%% doing filter image to have gradient magnitude
fx = imfilter(src_img, hx,'conv');
fy = imfilter(src_img, hy, 'conv');
% figure; imshow(fx );
% figure; imshow(fy);
mg = sqrt(fx.^2 + fy.^2);%magnitude of gradient vector

%% Non-maximum supression
%First, assuming that every pixels is non-edge
canny_edge_img = zeros(rows, cols);
phase_img = zeros(rows, cols);
% threshold = quantile(mg,0.5);
threshold = mean(mg(:));
count = 0;
for i = 2:rows-1
    for j = 2:cols-1
        if mg(i,j) > threshold       
            count = count + 1
            fxfy = fx(i,j) * fy(i,j);
            ep1 = 0;
            ep2 = 0;
            if fxfy == 0 %Case 1
                if fx == 0
                    ep1 = mg(i,j-1);
                    ep2 = mg(i,j+1);
                    phase_img(i,j) = 3;%90 degree
                elseif fy ==0
                    ep1 = mg(i-1,j);
                    ep2 = mg(i+1,j);   
                    phase_img(i,j) = 1;%0 degree
                end
            elseif(fxfy > 0)                
                if abs(fx) >= abs(fy) %Case 2
                    switch method
                        case 'interpolation'
                            alpha = abs(fy(i,j)/fx(i,j));
                            ep1 = mg(i+1,j+1) * alpha + mg(i+1,j) * (1-alpha);
                            ep2 = mg(i-1, j-1) * alpha + mg(i-1,j) * (1 - alpha);
                            phase_img(i,j) = 1;%0 degree
                        case 'average'
                            ep1 = (mg(i+1,j+1) + mg(i+1,j))/2;
                            ep2 = (mg(i-1, j-1) + mg(i-1,j))/2;
                            phase_img(i,j) = 1;%0 degree
                        case 'fast'
                            %Rouding to 0 degree
                            ep1 = mg(i-1,j);
                            ep2 = mg(i+1,j);
                            phase_img(i,j) = 1;%0 degree
                    end
                else %Case 3
                    switch method
                        case 'interpolation'
                            alpha = abs(fx(i,j)/fy(i,j));
                            ep1 = mg(i+1,j+1) * alpha + mg(i,j+1) * (1- alpha);
                            ep2 = mg(i-1,j-1) * alpha + mg(i-1,j) * (1 - alpha);
                            phase_img(i,j) = 2;%45 degree
                        case 'average'
                            ep1 = (mg(i+1,j+1) + mg(i,j+1))/2;
                            ep2 = (mg(i-1,j-1) + mg(i-1,j))/2;
                            phase_img(i,j) = 2;%45 degree
                        case 'fast'
                            %Rouding to 45 degree
                            ep1 = mg(i+1,j+1);
                            ep2 = mg(i-1,j-1);
                            phase_img(i,j) = 2;%45 degree
                    end
                    
                end                    
            elseif(fxfy < 0)                
                if abs(fx) >= abs(fy) %Case 4
                    switch method
                        case 'interpolation'
                            alpha = abs(fy(i,j)/fx(i,j));
                            ep1 = mg(i-1,j+1) * alpha + mg(i-1,j)*(1-alpha);
                            ep2 = mg(i+1,j-1) * alpha + mg(i+1,j) * (1-alpha);
                            phase_img(i,j) = 3;%90 degree
                        case 'average'
                            ep1 = (mg(i-1,j+1) + mg(i-1,j))/2;
                            ep2 = (mg(i+1,j-1) + mg(i+1,j))/2;
                            phase_img(i,j) = 3;%90 degree
                        case 'fast'
                            %Rouding to 90 degree
                            ep1 = mg(i,j+1);
                            ep2 = mg(i,j-1);
                            phase_img(i,j) = 3;%90 degree
                    end

                else %Case 5
                    switch method
                        case 'interpolation'
                            alpha = abs(fx(i,j)/fy(i,j));
                            ep1 = mg(i-1,j+1) * alpha + mg(i,j+1)*(1-alpha);
                            ep2 = mg(i+1,j-1) * alpha + mg(i,j-1) * (1-alpha);
                            phase_img(i,j) = 4;%135 degree
                        case 'average'
                            ep1 = (mg(i-1,j+1) + mg(i,j+1))/2;
                            ep2 = (mg(i+1,j-1) + mg(i,j-1))/2;
                            phase_img(i,j) = 4;%135 degree
                        case 'fast'
                            %Rouding to 135 degree
                            ep1 = mg(i-1,j+1);
                            ep2 = mg(i+1,j-1);
                            phase_img(i,j) = 4;%135 degree
                    end
                end                
            end

            if mg(i,j) >= max(ep1, ep2)
                %turn on
                canny_edge_img(i,j) = 1;
            end
        else
            %do nothing, it's not an edge point
        end
    end
end
phase_img = phase_img.*canny_edge_img; %remove under threshold points
%build color coded image
%yellow -> 0-45 degree
%green -> 45-90 degree
%blue -> 90-135 degree
%red -> 135-180 degree
color_edge_img = zeros(rows, cols, 3);
for x = 1:rows
    for y = 1:cols
        switch phase_img(x,y)
            case 1
                color_edge_img(x,y,1) = 255;
                color_edge_img(x,y,2) = 255;
            case 2
                color_edge_img(x,y,2) = 255;                
            case 3
                color_edge_img(x,y,3) = 255;
            case 4
                color_edge_img(x,y,1) = 255;
        end
    end
end
%% Post-processing
canny_edge_img = im2bw(canny_edge_img);
canny_edge_img = remove_isolated_point(canny_edge_img);
canny_edge_img = bwmorph(canny_edge_img,'thin');
end
